CHC LDA + MANOVA + centroids analysis
Dan Bruzzese
6/28/2020
using chc 6
present CV LDAs, have testing/training if need be
graphed LDAs are all the points.
MANOVA for species level
fit_full_species_man$pca_summary ## Importance of first k=7 (out of 35) components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.0413 0.7144 0.51001 0.38099 0.34918 0.30663 0.28705
## Proportion of Variance 0.3871 0.1822 0.09287 0.05183 0.04354 0.03357 0.02942
## Cumulative Proportion 0.3871 0.5694 0.66225 0.71408 0.75762 0.79119 0.82061summary(manova(as.matrix(data[,4:cols]) ~ hostRace * sex *site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 4 0.04507 53.399 28 1097.5 < 2.2e-16 ***
## sex 1 0.71370 17.421 7 304.0 < 2.2e-16 ***
## site 8 0.26950 8.085 56 1642.4 < 2.2e-16 ***
## hostRace:sex 4 0.80887 2.375 28 1097.5 8.147e-05 ***
## hostRace:site 1 0.96866 1.405 7 304.0 0.202640
## sex:site 7 0.77071 1.663 49 1547.8 0.002979 **
## hostRace:sex:site 1 0.97168 1.266 7 304.0 0.266990
## Residuals 310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PCA for Species Level
LDA for Species level
## Reference
## Prediction Cingulata Cornivora Mendax pom Zepheria
## Cingulata 29 0 0 0 0
## Cornivora 1 5 0 1 0
## Mendax 0 0 60 8 2
## pom 0 0 1 214 0
## Zepheria 0 0 0 1 15## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 9.584570e-01 9.213168e-01 9.312796e-01 9.771049e-01 6.646884e-01
## AccuracyPValue McnemarPValue
## 2.679678e-40 NaN
lda for species + sex + site
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 5.164179e-01 4.800429e-01 4.614714e-01 5.710717e-01 1.432836e-01
## AccuracyPValue McnemarPValue
## 4.176925e-58 NaN
Manova for pomonella
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * site * sex, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 1 0.70245 12.2841 7 203.00 4.504e-13 ***
## site 4 0.37641 8.1524 28 733.35 < 2.2e-16 ***
## sex 1 0.48494 30.8012 7 203.00 < 2.2e-16 ***
## hostRace:site 2 0.72182 5.1338 14 406.00 6.340e-09 ***
## hostRace:sex 1 0.95257 1.4440 7 203.00 0.189418
## site:sex 4 0.77940 1.8745 28 733.35 0.004284 **
## hostRace:site:sex 1 0.92648 2.3014 7 203.00 0.028089 *
## Residuals 209
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Pomonella
## Reference
## Prediction Apple Haw
## Apple 80 22
## Haw 26 96## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 7.857143e-01 5.693688e-01 7.261263e-01 8.375728e-01 5.267857e-01
## AccuracyPValue McnemarPValue
## 8.355647e-16 6.650055e-01
LDA predicting host and sex for pomonella
## Reference
## Prediction Apple_Female Apple_Male Haw_Female Haw_Male
## Apple_Female 51 5 17 2
## Apple_Male 1 21 1 8
## Haw_Female 19 1 44 5
## Haw_Male 0 8 4 37## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 6.830357e-01 5.662230e-01 6.177391e-01 7.434095e-01 3.169643e-01
## AccuracyPValue McnemarPValue
## 3.310270e-29 5.581411e-01
MANOVA for Zepheria
# no interaction because missing males at MtPleasant
summary(manova(as.matrix(data[,4:cols]) ~ sex + site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.74708 0.67708 5 10 0.6507
## site 1 0.83749 0.38808 5 10 0.8461
## Residuals 14LDA for Zepheria
## Reference
## Prediction Zepheria_Female Zepheria_Male
## Zepheria_Female 5 3
## Zepheria_Male 4 5## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.5882353 0.1793103 0.3292472 0.8155630 0.5294118
## AccuracyPValue McnemarPValue
## 0.4062810 1.0000000
LDA for sex +site zepheria
## Reference
## Prediction Zepheria_Female_EastLansing
## Zepheria_Female_EastLansing 2
## Zepheria_Female_MtPleasant 0
## Zepheria_Male_EastLansing 4
## Reference
## Prediction Zepheria_Female_MtPleasant
## Zepheria_Female_EastLansing 2
## Zepheria_Female_MtPleasant 1
## Zepheria_Male_EastLansing 0
## Reference
## Prediction Zepheria_Male_EastLansing
## Zepheria_Female_EastLansing 1
## Zepheria_Female_MtPleasant 5
## Zepheria_Male_EastLansing 2## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.29411765 -0.05699482 0.10313551 0.55958272 0.47058824
## AccuracyPValue McnemarPValue
## 0.95792652 0.03207164
MANOVA for Mendax
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.57727 6.1024 6 50 7.541e-05 ***
## site 2 0.10716 17.1232 12 100 < 2.2e-16 ***
## sex:site 2 0.74705 1.3082 12 100 0.2258
## Residuals 55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Mendax
## Reference
## Prediction Mendax_Female Mendax_Male
## Mendax_Female 20 15
## Mendax_Male 13 13## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.54098361 0.07072905 0.40849889 0.66935590 0.54098361
## AccuracyPValue McnemarPValue
## 0.55241652 0.85010674
LDA for sex +site mendax
## Reference
## Prediction Mendax_Female_Fenville Mendax_Female_OtisLake
## Mendax_Female_Fenville 7 1
## Mendax_Female_OtisLake 1 5
## Mendax_Female_Sewanee 0 2
## Mendax_Male_Fenville 6 0
## Mendax_Male_OtisLake 0 1
## Mendax_Male_Sewanee 0 0
## Reference
## Prediction Mendax_Female_Sewanee Mendax_Male_Fenville
## Mendax_Female_Fenville 0 5
## Mendax_Female_OtisLake 1 0
## Mendax_Female_Sewanee 4 1
## Mendax_Male_Fenville 0 5
## Mendax_Male_OtisLake 1 1
## Mendax_Male_Sewanee 4 0
## Reference
## Prediction Mendax_Male_OtisLake Mendax_Male_Sewanee
## Mendax_Female_Fenville 0 0
## Mendax_Female_OtisLake 1 2
## Mendax_Female_Sewanee 2 3
## Mendax_Male_Fenville 0 1
## Mendax_Male_OtisLake 1 1
## Mendax_Male_Sewanee 2 3## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.409836066 0.283523654 0.285504382 0.543223627 0.229508197
## AccuracyPValue McnemarPValue
## 0.001280557 NaN
MANOVA for cingulata
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.80383 1.2813 4 21 0.308955
## site 2 0.39036 3.1529 8 42 0.006927 **
## sex:site 2 0.66131 1.2059 8 42 0.319040
## Residuals 24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LDA for Cingulata
## Reference
## Prediction Cingulata_Female Cingulata_Male
## Cingulata_Female 8 6
## Cingulata_Male 7 9## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.5666667 0.1333333 0.3742735 0.7453925 0.5000000
## AccuracyPValue McnemarPValue
## 0.2923324 1.0000000
LDA for sex +site cingulata
## Reference
## Prediction Cingulata_Female_Fenville
## Cingulata_Female_Fenville 1
## Cingulata_Female_SouthBend 0
## Cingulata_Female_Urbana 0
## Cingulata_Male_Fenville 1
## Cingulata_Male_SouthBend 0
## Cingulata_Male_Urbana 1
## Reference
## Prediction Cingulata_Female_SouthBend Cingulata_Female_Urbana
## Cingulata_Female_Fenville 3 0
## Cingulata_Female_SouthBend 2 0
## Cingulata_Female_Urbana 4 1
## Cingulata_Male_Fenville 0 0
## Cingulata_Male_SouthBend 0 0
## Cingulata_Male_Urbana 1 1
## Reference
## Prediction Cingulata_Male_Fenville Cingulata_Male_SouthBend
## Cingulata_Female_Fenville 0 4
## Cingulata_Female_SouthBend 0 1
## Cingulata_Female_Urbana 0 1
## Cingulata_Male_Fenville 0 1
## Cingulata_Male_SouthBend 0 1
## Cingulata_Male_Urbana 0 2
## Reference
## Prediction Cingulata_Male_Urbana
## Cingulata_Female_Fenville 0
## Cingulata_Female_SouthBend 2
## Cingulata_Female_Urbana 0
## Cingulata_Male_Fenville 2
## Cingulata_Male_SouthBend 0
## Cingulata_Male_Urbana 0## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.17241379 0.04000000 0.05845608 0.35774755 0.34482759
## AccuracyPValue McnemarPValue
## 0.98826571 NaN
MANOVA for Cornivora
There is only one site
# no site because only sampled at one site.
summary(manova(as.matrix(data[,4:cols]) ~ sex, data = data), test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.54944 0.82004 2 2 0.5494
## Residuals 3LDA for cornivora
## Reference
## Prediction Cornivora_Female Cornivora_Male
## Cornivora_Female 0 2
## Cornivora_Male 0 2## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.50000000 0.00000000 0.06758599 0.93241401 1.00000000
## AccuracyPValue McnemarPValue
## 1.00000000 0.47950012
lines of diversification
euclidian distance
## Apple_Female Apple_Male Cingulata_Female Cingulata_Male
## Apple_Female 0.0000000
## Apple_Male 1.1052241 0.0000000
## Cingulata_Female 1.4738114 1.6036342 0.0000000
## Cingulata_Male 1.4051571 1.6128575 0.1130547 0.0000000
## Cornivora_Female 1.8164505 1.6514454 0.4825512 0.5956054
## Cornivora_Male 2.3568964 1.7473686 1.3312107 1.4407462
## Haw_Female 0.5292825 1.5234184 1.9508336 1.8698619
## Haw_Male 0.6482415 0.5290546 1.6459586 1.6184269
## Mendax_Female 1.1742028 2.0634150 1.2374012 1.1243501
## Mendax_Male 1.0269548 1.8389388 1.0200228 0.9078379
## Zepheria_Female 2.2810030 3.3164335 2.4854319 2.3757562
## Zepheria_Male 2.3942415 3.4468577 2.6586999 2.5484713
## Cornivora_Female Cornivora_Male Haw_Female Haw_Male
## Apple_Female
## Apple_Male
## Cingulata_Female
## Cingulata_Male
## Cornivora_Female 0.0000000
## Cornivora_Male 0.8731903 0.0000000
## Haw_Female 2.3257992 2.8860880 0.0000000
## Haw_Male 1.8371458 2.1364704 1.0022950 0.0000000
## Mendax_Female 1.7199319 2.5456134 1.3364646 1.7570533
## Mendax_Male 1.5002006 2.3102051 1.2814072 1.5661386
## Zepheria_Female 2.9567427 3.8164525 2.1620829 2.9257376
## Zepheria_Male 3.1317505 3.9891585 2.2408291 3.0416391
## Mendax_Female Mendax_Male Zepheria_Female Zepheria_Male
## Apple_Female
## Apple_Male
## Cingulata_Female
## Cingulata_Male
## Cornivora_Female
## Cornivora_Male
## Haw_Female
## Haw_Male
## Mendax_Female 0.0000000
## Mendax_Male 0.2497539 0.0000000
## Zepheria_Female 1.3207469 1.5704754 0.0000000
## Zepheria_Male 1.4794064 1.7289277 0.1809107 0.0000000mds plot
centroids
centroid slopes
** should I calc slope for each host race pair? or only two groups as shown below?
pomonella == y ~ -0.8633 + 1.4937 * LD2.mean
rest == y ~ 1.8431 + 1.5561 * LD2.mean
these slopes are very similar. Just different intercepts. relationship btwn axes is the same for all species, but the pomonella host races have greater LD2 values compared to the other species. Also female to male differentiation flows along the same slope except the relationship is reversed in zepheria and cingulata. Cingulata, zepheria, and mendax have small sex diff. Pomonella has great sex differences. Cornivora look to have great sex differences but further sampling will be needed to confrim.
##
## Call:
## lm(formula = LD1.mean ~ LD2.mean, data = .)
##
## Residuals:
## 1 2 3 4
## 0.40664 -0.10011 -0.08027 -0.22626
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.8633 0.2246 -3.843 0.0615 .
## LD2.mean 1.4937 0.5724 2.609 0.1208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3413 on 2 degrees of freedom
## Multiple R-squared: 0.7729, Adjusted R-squared: 0.6594
## F-statistic: 6.808 on 1 and 2 DF, p-value: 0.1208##
## Call:
## lm(formula = LD1.mean ~ LD2.mean, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48301 -0.15745 0.06621 0.23799 0.31705
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8431 0.1729 10.659 4.02e-05 ***
## LD2.mean 1.5561 0.1606 9.686 6.94e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3161 on 6 degrees of freedom
## Multiple R-squared: 0.9399, Adjusted R-squared: 0.9299
## F-statistic: 93.83 on 1 and 6 DF, p-value: 6.944e-05mcmc sampling ???
calc angle
theta = tan_inverse (slope)
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)
I think the slope is given below in radians? y1, x1= always the female and y0 an x0 are always males.it should’t matter though…
# apple
ang = centroids %>% filter(host_Species == "Apple")
ang## host_Species sex LD1.mean LD2.mean
## 1 Apple Female -0.115682 0.2282819
## 2 Apple Male -1.160544 -0.1319735y1 = ang$LD2.mean[1] # female y
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x
x0 = ang$LD1.mean[2] #male x
slope = (y1 - y0)/(x1 - x0)
slope## [1] 0.3447876theta = atan(slope)
theta## [1] 0.3320237apple = cbind(slope,theta)
apple = as.data.frame(apple)
apple$host = "apple"
apple## slope theta host
## 1 0.3447876 0.3320237 apple# haw
ang = centroids %>% filter(host_Species == "Haw")
ang## host_Species sex LD1.mean LD2.mean
## 1 Haw Female 0.1114762 0.7063395
## 2 Haw Male -0.7638430 0.2180654y1 = ang$LD2.mean[1] # female y
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x
x0 = ang$LD1.mean[2] #male x
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)
haw = cbind(slope,theta)
haw = as.data.frame(haw)
haw$host = "haw"
haw## slope theta host
## 1 0.5578241 0.5088304 haw# zepheria
ang = centroids %>% filter(host_Species == "Zepheria")
ang## host_Species sex LD1.mean LD2.mean
## 1 Zepheria Female 2.153432 -0.004308698
## 2 Zepheria Male 2.276472 0.128317729y1 = ang$LD2.mean[1] # female y
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x
x0 = ang$LD1.mean[2] #male x
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)
Zepheria = cbind(slope,theta)
Zepheria = as.data.frame(Zepheria)
Zepheria$host = "Zepheria"
Zepheria## slope theta host
## 1 1.077911 0.8228755 Zepheria# mendax
ang = centroids %>% filter(host_Species == "Mendax")
ang## host_Species sex LD1.mean LD2.mean
## 1 Mendax Female 0.8877174 -0.3815885
## 2 Mendax Male 0.6495075 -0.4566423y1 = ang$LD2.mean[1] # female y
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x
x0 = ang$LD1.mean[2] #male x
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)
Mendax = cbind(slope,theta)
Mendax = as.data.frame(Mendax)
Mendax$host = "Mendax"
Mendax## slope theta host
## 1 0.3150745 0.3052286 Mendax# cingulata
ang = centroids %>% filter(host_Species == "Cingulata")
ang## host_Species sex LD1.mean LD2.mean
## 1 Cingulata Female -0.002386062 -1.241168
## 2 Cingulata Male 0.079528433 -1.163249y1 = ang$LD2.mean[1] # female y
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x
x0 = ang$LD1.mean[2] #male x
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)
Cingulata = cbind(slope,theta)
Cingulata = as.data.frame(Cingulata)
Cingulata$host = "Cingulata"
Cingulata## slope theta host
## 1 0.9512238 0.7604056 Cingulata# Cornivora
ang = centroids %>% filter(host_Species == "Cornivora")
ang## host_Species sex LD1.mean LD2.mean
## 1 Cornivora Female -0.353129 -1.572582
## 2 Cornivora Male -1.170673 -1.879313y1 = ang$LD2.mean[1] # female y
y0 = ang$LD2.mean[2] #male y
x1 = ang$LD1.mean[1] # female x
x0 = ang$LD1.mean[2] #male x
slope = (y1 - y0)/(x1 - x0)
theta = atan(slope)
Cornivora = cbind(slope,theta)
Cornivora = as.data.frame(Cornivora)
Cornivora$host = "Cornivora"
Cornivora## slope theta host
## 1 0.3751856 0.3589334 Cornivorarbind(apple, haw, Zepheria, Mendax, Cingulata, Cornivora)## slope theta host
## 1 0.3447876 0.3320237 apple
## 2 0.5578241 0.5088304 haw
## 3 1.0779109 0.8228755 Zepheria
## 4 0.3150745 0.3052286 Mendax
## 5 0.9512238 0.7604056 Cingulata
## 6 0.3751856 0.3589334 Cornivorasite centroids
